.libPaths(c("/data/home/lyang/R/x86_64-redhat-linux-gnu-library/4.1",
            "/data/home/bioinfo/R/R4.1.0/library_B3.13",
            "/usr/lib64/R/library","/usr/share/R/library") )
dyn.load("/data/home/bioinfo/programs/hdf5-1.12.0/hdf5/lib/libhdf5_hl.so.200")
knitr::opts_knit$set(root.dir = "/data/home/lyang/Visium_spotlight")
setwd("/data/home/lyang/Visium_spotlight")

Load packages

# library(Matrix)
library(data.table)
library(Seurat)
library(SeuratDisk)
library(dplyr)


dataset_name = "cell2location34"
# in_house  cell2location34
sc_dataset_file = paste(dataset_name,"sc_dataset.rds",sep ="_")
rds_file= sc_dataset_file
sc_seu <- readRDS(rds_file)
  
rds_file='Visium_spatial.rds'
sp_seu <- readRDS(rds_file)

Quick data exploration, cell number of each cell type:

# "cell2location34" Tabula in_house
dataset_name = "cell2location34"
sc_dataset = sc_seu
if(dataset_name == "cell2location34"){

  cell_type_table = as.data.frame(table(sc_dataset$Subset))

}else if(dataset_name == "in_house"){
  cell_type_table = as.data.frame(table(sc_dataset$blue.main))
}else if(dataset_name == "cell2location44"){
  cell_type_table = as.data.frame(table(sc_dataset$CellType))
}




colnames(cell_type_table) = c("cell_type","frequency")
library(ggplot2)
# Basic barplot
p <- ggplot(data=cell_type_table, aes(x=cell_type, y=frequency,fill=cell_type)) +
  geom_bar(stat="identity")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))+
       geom_text(aes(label=frequency), position=position_dodge(width=0.9), vjust=-0.25)


# ggsave(p,filename = "cell_type_histogram.pdf",width = 15,
#        height = 20)

print(p)

getwd()
## [1] "/data/home/lyang/Visium_spotlight"
if(file.exists("old_version_ls.rds")){
  rds_file= "old_version_ls.rds"
  spotlight_ls <- readRDS(rds_file)
  nmf_mod <- spotlight_ls[[1]]
  decon_mtrx <- spotlight_ls[[2]]
  }else{
    Idents(sc_seu) = sc_seu$Subset
sc_seu <- SCTransform(sc_seu, verbose = FALSE) %>% RunPCA(verbose = FALSE)
sp_seu <- SCTransform(sp_seu,assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)

#### Extract the top marker genes from each cluster ####
Seurat::Idents(object = sc_seu) <- sc_seu@meta.data$Subset
cluster_markers_all <- Seurat::FindAllMarkers(object = sc_seu, 
                                              assay = "RNA",
                                              slot = "data",
                                              verbose = TRUE, 
                                              only.pos = TRUE, 
                                              logfc.threshold = 1,
                                              min.pct = 0.9)

saveRDS(cluster_markers_all,"cluster_markers_all.rds")
setwd("/data/home/lyang/Visium_spotlight")


# count_matrix <- as.matrix(sc_seu@assays$RNA@counts)
# 
# # removing genes not expressed in 3% of the cells
# min_cells <- 0.03 * ncol(sc_seu)
# 
# logic_vector = rowSums(count_matrix != 0) > min_cells
# 
# abundant_genes = rownames(sc_seu@assays$RNA@counts[logic_vector,])
# sc_seu <- subset(sc_seu, features = abundant_genes)

set.seed(123)
start_time <- Sys.time()
library(SPOTlight, lib.loc="/data/home/bioinfo/R/R4.1.0/library_B3.13/old_versions/")
spotlight_ls <- spotlight_deconvolution(se_sc = sc_seu,
                                      counts_spatial = sp_seu@assays$Spatial@counts,
                                      clust_vr = "Subset",
                                      cluster_markers = cluster_markers_all,
                                      cl_n = 100, # 100 by default
                                      hvg = 5000,
                                      ntop = NULL,
                                      transf = "uv",
                                      method = "nsNMF",
                                      min_cont = 0.09)
saveRDS(spotlight_ls,"old_version_ls.rds")
end_time <- Sys.time()
end_time - start_time
    
  }
unloadNamespace("SPOTlight")
library(SPOTlight, lib.loc="/data/home/bioinfo/R/R4.1.0/library_B3.13/old_versions/")
h <- NMF::coef(nmf_mod[[1]])
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
  h = h,
  train_cell_clust = nmf_mod[[2]])

topic_profile_plts[[2]] + ggplot2::theme(
  axis.text.x = ggplot2::element_text(angle = 90), 
  axis.text = ggplot2::element_text(size = 12))

topic_profile_plts[[1]] + theme(axis.text.x = element_text(angle = 90), 
                                axis.text = element_text(size = 12))

# change the cell type names in sp_seu, replace _ with . to keep names coincide between sp_seu and decon_mtrx
decon_mtrx_sub <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
# decon_mtrx_sub[decon_mtrx_sub < 0.08] <- 0
decon_mtrx <- cbind(decon_mtrx_sub, "res_ss" = decon_mtrx[, "res_ss"])
rownames(decon_mtrx) <- colnames(sp_seu)

decon_df <- decon_mtrx %>%
  data.frame() %>%
  tibble::rownames_to_column("barcodes")

sp_seu@meta.data <- sp_seu@meta.data %>%
  tibble::rownames_to_column("barcodes") %>%
  dplyr::left_join(decon_df, by = "barcodes") %>%
  tibble::column_to_rownames("barcodes")

norm_weights_matrix = decon_mtrx

# Seurat::SpatialFeaturePlot(
#   object = sp_seu,
#   features = c("B.Cycling", "B.GC.DZ"),
#   alpha = c(0.1, 1))
cell_types_all <- colnames(decon_mtrx)[which(colnames(decon_mtrx) != "res_ss")]
getwd()
## [1] "/data/home/lyang/Visium_spotlight"
# install.packages("imager")
# SPOTlight::spatial_scatterpie(se_obj = sp_seu,
#                               cell_types_all = cell_types_all,
#                               img_path = "tissue_lowres_image.png",
#                               pie_scale = 0.4)
# 
# 
# SPOTlight::spatial_scatterpie(se_obj = sp_seu,
#                               cell_types_all = cell_types_all,
#                               img_path = "tissue_lowres_image.png",
#                               cell_types_interest = "B.naive",
#                               pie_scale = 0.4)
# "T_CD4+_naive","FDC","B_naive"
# This can be dynamically visualized with DT as shown below
unloadNamespace("SPOTlight")
library(SPOTlight)

ct <- colnames(decon_mtrx[,-ncol(decon_mtrx)])
# decon_mtrx[decon_mtrx < 0.1] <- 0


library(RColorBrewer)

if(dataset_name == "Tabula"){
  n <- 29
}else if(dataset_name == "cell2location"){
  n <- 44
}else if(dataset_name == "cell2location34"){
  n <- 34
}

  
colrs <- brewer.pal.info[brewer.pal.info$colorblind == TRUE, ]
col_vec = unlist(mapply(brewer.pal, colrs$maxcolors, rownames(colrs)))
set.seed(123)
col <- sample(col_vec, n)


# show T cell and 
# B cell zones and GCs with follicular dendritic cells (FDCs)
# FDC, T_CD4+_naive, B_naive
# Transcriptionally fine subtypes (B_GC_DZ, B_GC_LZ, B_GC_prePB and 
# T_CD4+_TfH_GC); transcriptionally distinct subtypes (B_Cycling and FDC)


paletteMartin <- col

pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
pal_back <- pal


plot_3_region <- function(pal)
{
  for (char in names(pal)) {
  # print(char)
  if(char %in% c("T.CD4..naive","FDC","B.naive") ){
    if(char == "T.CD4..naive")
    {
      pal[char] = "#FFFF00"
    } else if (char == "FDC")
    {
      pal[char] = "#add8e6"
    } else if (char == "B.naive")
    {
      pal[char] = "#FF0000"
    }
    next
  }
  pal[char] = "#00008b"
}

  return(pal)
}
packageVersion("SPOTlight")
## [1] '0.99.4'
pal = plot_3_region(pal_back)
library(ggplot2)
plotSpatialScatterpie(
  x = sp_seu,
  y = decon_mtrx,
  cell_types = colnames(y),
  img = FALSE,
  scatterpie_alpha = 1,
  pie_scale = 0.4) +
  scale_fill_manual(
    values = pal,
    breaks = names(pal))+ scale_y_reverse()

plot_n_cell_type <- function(pal,col_vec,n, cell_type_vector)
{
  for (char in names(pal)) {
  # print(char)
  if(char %in%  cell_type_vector){
    next
  }
  pal[char] = "#00008b"
}

  return(pal)
}

cell_type_vector = c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC",
                 "B.Cycling",  "FDC")

pal = plot_n_cell_type(pal_back,col_vec,6,cell_type_vector)
plotSpatialScatterpie(
  x = sp_seu,
  y = decon_mtrx,
  cell_types = colnames(y),
  img = FALSE,
  scatterpie_alpha = 1,
  pie_scale = 0.4) +
  scale_fill_manual(
    values = pal,
    breaks = names(pal))+ scale_y_reverse()

pal = pal_back
# type(pal)
plotSpatialScatterpie(
  x = sp_seu,
  y = decon_mtrx,
  cell_types = colnames(y),
  img = FALSE,
  scatterpie_alpha = 1,
  pie_scale = 0.4) +
  scale_fill_manual(
    values = pal,
    breaks = names(pal))+ scale_y_reverse()

rds_file='predictions.assay.rds'
predictions.assay <- readRDS(rds_file)
mat = decon_mtrx
spatial = sp_seu
predictions.assay@data <- t(mat)
meta.features <- as.data.frame(colnames(mat) )
rownames(meta.features) = meta.features[,1]
meta.features$`colnames(mat)` =NULL
predictions.assay@meta.features = meta.features


spatial[["predictions"]] <- predictions.assay

DefaultAssay(spatial) <- "predictions"
# table(sc_dataset@meta.data$Subset)


if(dataset_name == "Tabula"){
  
p = SpatialFeaturePlot(spatial, features = c("stromal cell", "t cell"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# print(p)
# p = SpatialFeaturePlot(spatial, features = c("naive b cell", "b cell"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# print(p)
  
}else if(dataset_name == "cell2location"){
  
 SpatialFeaturePlot(spatial, features = c("T_Treg", "DC_cDC2"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# SpatialFeaturePlot(spatial, features = c("DC_cDC1", "DC_cDC2"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)

}else if(dataset_name == "cell2location34"){
  
# detach("package:SpatialExperiment", unload=TRUE)
  
p = SpatialFeaturePlot(spatial, features = c("B.GC.LZ", "T.CD4..TfH.GC","B.GC.prePB","FDC"), pt.size.factor = 1.6, ncol = 4, crop = TRUE) + ggtitle("germinal center light zone")
print(p)
 # + scale_fill_continuous(limits = c(0, 1))

p = SpatialFeaturePlot(spatial, features = c("B.Cycling", "B.GC.DZ"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)+ ggtitle("germinal center dark zone") 
print(p)
SpatialFeaturePlot(spatial, features = c("B.naive", "B.preGC"), pt.size.factor = 1.6, ncol = 2, crop = TRUE) + ggtitle("B follicle + pre GC") 
}

# table(sc_dataset@meta.data$Subset)
# https://github.com/satijalab/seurat/blob/master/R/visualization.R
# Leave out selected cell types from sc reference dataset(B_GC_DZ, B_GC_LZ, B_GC_prePB, T_CD4+_TfH_GC, B_Cycling and FDC, these six subtypes are expected to be present in annotated GC (positive locations) one by one using ‘subset’ command.


# sc_seu <- SCTransform(sc_seu, verbose = FALSE) %>% RunPCA(verbose = FALSE)
# sp_seu <- SCTransform(sp_seu,assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)
# saveRDS(sp_seu,"sp_seu_SCTransform.rds")
# saveRDS(sc_seu,"sc_seu_SCTransform.rds")


rds_file= "sp_seu_SCTransform.rds"
sp_seu <- readRDS(rds_file)
rds_file= "sc_seu_SCTransform.rds"
sc_seu <- readRDS(rds_file)

#### Extract the top marker genes from each cluster ####
# Seurat::Idents(object = sc_seu) <- sc_seu@meta.data$Subset
# cluster_markers_all <- Seurat::FindAllMarkers(object = sc_seu, 
#                                               assay = "RNA",
#                                               slot = "data",
#                                               verbose = TRUE, 
#                                               only.pos = TRUE, 
#                                               logfc.threshold = 1,
#                                               min.pct = 0.9)
# saveRDS(cluster_markers_all,"cluster_markers_all.rds")
# Leave out selected cell types from  marker genes variables

rds_file= "cluster_markers_all.rds"
cluster_markers_all <- readRDS(rds_file)

#  unloadNamespace("SPOTlight") 
# library(SPOTlight, lib.loc="/data/home/bioinfo/R/R4.1.0/library_B3.13/old_versions/")
# for (selected_type in c( "B_GC_LZ", "B_GC_prePB","T_CD4+_TfH_GC",
#                  "B_Cycling",  "FDC") )
# {
#   # selected_type = "B_GC_DZ"
# 
#   spotlight_result_file = paste(dataset_name,"spotlight_deleted", selected_type, ".rds",sep ="_")
#   sc_reference <- subset(sc_seu, subset = Subset != selected_type)
#   
#   #### Extract the top marker genes from each cluster ####
#   Seurat::Idents(object = sc_reference) <- sc_reference@meta.data$Subset
# 
#   set.seed(123)
#   setwd("/data/home/lyang/Visium_spotlight")
# 
#   spotlight_ls <- spotlight_deconvolution(se_sc = sc_reference,
#                                         counts_spatial = sp_seu@assays$Spatial@counts,
#                                         clust_vr = "Subset",
#                                         cluster_markers = cluster_markers_all,
#                                         cl_n = 100, # 100 by default
#                                         hvg = 5000,
#                                         ntop = NULL,
#                                         transf = "uv",
#                                         method = "nsNMF",
#                                         min_cont = 0.09)
# 
#   saveRDS(spotlight_ls,spotlight_result_file)
#   print(spotlight_result_file)
#     
# }
# getwd()
rmsd_summary = rep(0,6)
names(rmsd_summary) =  c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC",
                 "B.Cycling",  "FDC")

manual_GC_annotation = read.csv("/data/home/lyang/scripts_git/manual_GC_annotation.csv")
rds_file = "Visium_spatial.rds"
spatial <- readRDS(rds_file)
spatial@meta.data$Barcode = rownames(spatial@meta.data)
spatial@meta.data = merge(spatial@meta.data, manual_GC_annotation, by = "Barcode", sort = FALSE)

spatial@meta.data[spatial@meta.data$GC == '',"GC"] = "no"

p2 <- SpatialDimPlot(spatial, group.by = "GC") 
# p1 + p2
print(p2)

for (selected_type in  c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC",
                 "B.Cycling",  "FDC") )
{
  # selected_type = "B.GC.DZ"

  spotlight_result_file = paste(dataset_name,"spotlight_deleted", selected_type, ".rds",sep ="_")
  setwd("/data/home/lyang/Visium_spotlight")
  rds_file= spotlight_result_file
  spotlight_ls <- readRDS(rds_file)
  nmf_mod <- spotlight_ls[[1]]
  decon_mtrx <- spotlight_ls[[2]]
  
  
  
  # change the cell type names in sp_seu, replace _ with . to keep names coincide between sp_seu and decon_mtrx
  decon_mtrx_sub <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
  # decon_mtrx_sub[decon_mtrx_sub < 0.08] <- 0
  decon_mtrx <- cbind(decon_mtrx_sub, "res_ss" = decon_mtrx[, "res_ss"])
  rownames(decon_mtrx) <- colnames(sp_seu)
  
  decon_df <- decon_mtrx %>%
    data.frame() %>%
    tibble::rownames_to_column("barcodes")
  
  sp_seu@meta.data <- sp_seu@meta.data %>%
    tibble::rownames_to_column("barcodes") %>%
    dplyr::left_join(decon_df, by = "barcodes") %>%
    tibble::column_to_rownames("barcodes")  
  # delete last column that is red_ss value
  matrix_result_complete = norm_weights_matrix[,-ncol(norm_weights_matrix)]

  matrix_result_deleted = as.data.frame(decon_mtrx[,-ncol(decon_mtrx)])
  # matrix_result_deleted lack one column that is the deleted cell type's column,
  # we add this deleted cell type's column after last cloumn
  matrix_result_deleted$selected_type = rep(0,nrow(matrix_result_deleted))
  #  meanwhile we cut the  deleted cell type's column in complete matrix_result_complete and insert it into matrix_result_complete's last column so deleted matrix and complete matrix have same cell type column oder 
  selected_column = matrix_result_complete[,selected_type]
  matrix_result_complete = as.data.frame(matrix_result_complete)
  df = matrix_result_complete[,!(names(matrix_result_complete) == selected_type)]
  df$selected_type = selected_column
  matrix_result_complete = as.matrix(df)
  matrix_result_deleted = as.matrix(matrix_result_deleted)
  
  
  # for each spatial spot, Calculate the root-mean-square error (RMSE) distance between the predicted proportion under different reference dataset situations
  # load("/data/home/lyang/Visium_RCTD/.RData")
  
  
  matrix_substrated = matrix_result_complete - matrix_result_deleted
  rmsd_per_spot = sqrt(rowSums(matrix_substrated^2)/ncol(matrix_substrated))
  rmsd_per_situation = sum(rmsd_per_spot)
  rmsd_summary[selected_type] = rmsd_per_situation
  
  
  
  
  
  
  
  
  
  unloadNamespace("SPOTlight")
  library(SPOTlight)
  
  ct <- colnames(decon_mtrx)
  # decon_mtrx[decon_mtrx < 0.1] <- 0
  
  
  # # show T cell and 
  # # B cell zones and GCs with follicular dendritic cells (FDCs)
  # # FDC, T_CD4+_naive, B_naive
  # # Transcriptionally fine subtypes (B_GC_DZ, B_GC_LZ, B_GC_prePB and 
  # # T_CD4+_TfH_GC); transcriptionally distinct subtypes (B_Cycling and FDC)

  library(ggplot2)
  
  plot_n_cell_type <- function(pal,col_vec,n, cell_type_vector)
  {
    for (char in names(pal)) {
    # print(char)
    if(char %in%  cell_type_vector){
      next
    }
    pal[char] = "#00008b"
  }
  
    return(pal)
  }
  
  cell_type_vector = c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC",
                   "B.Cycling",  "FDC")
  
  pal = plot_n_cell_type(pal_back,col_vec,6,cell_type_vector)
  p1 = plotSpatialScatterpie(
    x = sp_seu,
    y = decon_mtrx[,-ncol(decon_mtrx)],
    cell_types = colnames(y),
    img = FALSE,
    scatterpie_alpha = 1,
    pie_scale = 0.4) +
    scale_fill_manual(
      values = pal,
      breaks = names(pal)) + scale_y_reverse()+ggtitle(selected_type)
   print(p1)
}

df <- data.frame(deleted_cell_type=names(rmsd_summary),
                RMSD=rmsd_summary  )
p<-ggplot(data=df, aes(x=deleted_cell_type, y=RMSD)) +
  geom_bar(stat="identity")
p

# Evaluate the area under Precision Recall curve separately for each cell type, then treat all classes equally and take the average across remaining 5 GC cell types
library(ROCR)
average_PR_curve.area = rep(0,7)
names(average_PR_curve.area) = c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC", "B.Cycling",  "FDC","complete")
manual_GC_annotation = read.csv("/data/home/lyang/scripts_git/manual_GC_annotation.csv")

getwd()
## [1] "/data/home/lyang/Visium_spotlight"
for (deleted_type in c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC", "B.Cycling",  "FDC","complete") )
{

  if(deleted_type == "complete"){
    rds_file= "old_version_ls.rds"
  }else{
  spotlight_result_file = paste(dataset_name,"spotlight_deleted", deleted_type, ".rds",sep ="_")
  rds_file= spotlight_result_file
  }
    spotlight_ls <- readRDS(rds_file)
    nmf_mod <- spotlight_ls[[1]]
    decon_mtrx <- spotlight_ls[[2]]  
    
    
  # change the cell type names in sp_seu, replace _ with . to keep names coincide between sp_seu and decon_mtrx
  decon_mtrx_sub <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
  # decon_mtrx_sub[decon_mtrx_sub < 0.08] <- 0
  decon_mtrx <- cbind(decon_mtrx_sub, "res_ss" = decon_mtrx[, "res_ss"])
  rownames(decon_mtrx) <- colnames(sp_seu)
  
  decon_df <- decon_mtrx %>%
    data.frame() %>%
    tibble::rownames_to_column("barcodes")
  
  sp_seu@meta.data <- sp_seu@meta.data %>%
    tibble::rownames_to_column("barcodes") %>%
    dplyr::left_join(decon_df, by = "barcodes") %>%
    tibble::column_to_rownames("barcodes")  

  matrix_result_deleted = decon_mtrx[,-ncol(decon_mtrx)]
  sum_PR_curve.area = 0
  # Evaluate the area under Precision Recall curve separately for each cell type
  i = 1
  for (remain_type in c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC", "B.Cycling",  "FDC") )
{
    if(remain_type == deleted_type){
      next
    }
      # print(c("deleted_type",deleted_type,"remain_type",remain_type))
      # extract remain_type's predicted proportions in each spot
    proportions = as.data.frame(matrix_result_deleted[,remain_type]) 
    colnames(proportions) = "proportion"
    # set GC spot's value to 1, non-GC spot's value to 0
    manual_GC_annotation[manual_GC_annotation$GC == '',"GC"] = 0
    manual_GC_annotation[manual_GC_annotation$GC == 'GC',"GC"] = 1
    proportions$Barcode = rownames(proportions)
    proportions_label = merge(proportions,manual_GC_annotation,by = "Barcode")
 
    library(ggplot2)   

    if(remain_type != "FDC"){
        print(paste("deleted_type:",deleted_type,"remain_type:",remain_type,
            "predicted_positive:",sum(proportions_label$proportion>0),
            "actual_positive:",sum(proportions_label$GC>0)))
      
        print("GC spot")
        print(table(proportions_label[proportions_label$GC==1,"proportion"]))
        print("non-GC spot")
        print(table(proportions_label[proportions_label$GC==0,"proportion"]))
    p3 = ggplot(proportions_label[proportions_label$proportion>0,], aes(x=proportion))+
    geom_histogram(color="darkblue", fill="lightblue") + facet_grid(. ~ GC)+
      ggtitle(paste("deleted_type:",deleted_type,"remain_type:",remain_type,
            "predicted_positive:",sum(proportions_label$proportion>0),
            "actual_positive:",sum(proportions_label$GC>0)))
    print(p3)

    }
    
    
    
    pred <- prediction(proportions_label[,2],proportions_label[,3])
    
    PR_curve.area = performance(pred, measure = "aucpr")
    PR_curve.area@y.values[[1]]
    sum_PR_curve.area = PR_curve.area@y.values[[1]] + sum_PR_curve.area
    
    PR_curve = performance(pred, measure = "prec", x.measure = "rec")
    if(i ==1)
  {
  df_all <- data.frame(cell_type=rep(remain_type,each=length(PR_curve@x.values[[1]])), 
  recall=c(PR_curve@x.values[[1]]),
  precision=c(PR_curve@y.values[[1]]),
  cutoff = c(PR_curve@alpha.values[[1]])) 
  i = i + 1 
    }else{

  df <- data.frame(cell_type=rep(remain_type,each=length(PR_curve@x.values[[1]])), 
  recall=c(PR_curve@x.values[[1]]),
  precision=c(PR_curve@y.values[[1]]),
  cutoff = c(PR_curve@alpha.values[[1]]))
  df_all =  rbind(df_all,df)    
  i = i + 1 
  }
  }

  df_all$cell_type = as.factor(df_all$cell_type)
  plt <- ggplot(df_all, aes(x=recall, y=precision, color=cell_type))+geom_point() +
    ggtitle(deleted_type)
  print(plt)
  
  average_PR_curve.area[deleted_type] = sum_PR_curve.area/5

  
}
## [1] "deleted_type: B.GC.DZ remain_type: B.GC.LZ predicted_positive: 49 actual_positive: 378"
## [1] "GC spot"
## 
##                 0   0.1577099103182 0.163724812431107 0.168858996150598 
##               332                 1                 1                 1 
## 0.172351245858612 0.173193106770393 0.185556643998863 0.187939123313963 
##                 1                 1                 1                 1 
## 0.190502501322252 0.194593680081007 0.196856265886727 0.199745974191692 
##                 1                 1                 1                 1 
## 0.202901818496328 0.205375191473186 0.206871892833241 0.208639504533753 
##                 1                 1                 1                 1 
## 0.209536971974085 0.213379008205849 0.216359328575328 0.217588628923261 
##                 1                 1                 1                 1 
## 0.222190585429564 0.223517787889275 0.223622915803676 0.226098595480703 
##                 1                 1                 1                 1 
## 0.228080234025031 0.230106203619952 0.233569095706999 0.233982259409767 
##                 1                 1                 1                 1 
## 0.239617756805643 0.240553575183547 0.240749281335039 0.241250305972474 
##                 1                 1                 1                 1 
## 0.241749733383228 0.244203269424668 0.251745063287218 0.253676284730386 
##                 1                 1                 1                 1 
## 0.261401099367808   0.2623247256102 0.263561129502264 0.266618413992399 
##                 1                 1                 1                 1 
##  0.27407628390125 0.278347191556542 0.287577296269503 0.292408599653626 
##                 1                 1                 1                 1 
## 0.305389592270722 0.307488903934843  0.32825944300301 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.150323761753507  0.18087053777631 0.267945587977937 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.DZ remain_type: B.GC.prePB predicted_positive: 3 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.184111852907085 
##               377                 1 
## [1] "non-GC spot"
## 
##                 0 0.119256501468096 0.146838209587556 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.DZ remain_type: T.CD4..TfH.GC predicted_positive: 7 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.181125588997005 0.232981432261497 0.250090970364881 
##               375                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.121445880897914 0.161786132870309 0.210907105936577 
##              3653                 1                 1                 1 
## 0.316738167588768 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.DZ remain_type: B.Cycling predicted_positive: 40 actual_positive: 378"
## [1] "GC spot"
## 
##                 0  0.14775060078779 0.152815147705166  0.16164777460797 
##               339                 1                 1                 1 
## 0.165984386962632 0.171350456299436 0.179602840543764 0.181036287951638 
##                 1                 1                 1                 1 
## 0.188007698502191 0.197706411823626 0.198807261625983 0.202156476398911 
##                 1                 1                 1                 1 
## 0.202254303156224 0.204803117328647 0.205375257529594 0.206317562071278 
##                 1                 1                 1                 1 
## 0.206371603305384 0.207155774601566 0.212003506449966 0.213954108074288 
##                 1                 1                 1                 1 
## 0.216223087206259 0.221346670159729 0.223167094866398 0.225083919381899 
##                 1                 1                 1                 1 
## 0.225653926522961 0.226048649231645 0.233124045526933 0.234015885860872 
##                 1                 1                 1                 1 
## 0.236824163627549  0.24017024198831 0.244789460685742 0.244873105907738 
##                 1                 1                 1                 1 
## 0.260409060312816 0.263487790509182 0.263957884044737 0.274921247507442 
##                 1                 1                 1                 1 
## 0.288024738299017 0.288782645670002 0.346772744405923 0.438167560902959 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.187030656023299 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: B.GC.LZ remain_type: B.GC.DZ predicted_positive: 7 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.157455394398931 0.187712539963545 0.199863323019129 
##               373                 1                 1                 1 
## 0.214085202414908 0.234853696262563 
##                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.172644309401411 0.251507496943164 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.LZ remain_type: B.GC.prePB predicted_positive: 5 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.229566249898306 0.234127932368522 
##               376                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.127379210894173  0.14098119183902 0.260707072615927 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.LZ remain_type: T.CD4..TfH.GC predicted_positive: 4 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.173767313023636 
##               377                 1 
## [1] "non-GC spot"
## 
##                 0 0.232768492251587 0.236930610648794 0.463681078063691 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.LZ remain_type: B.Cycling predicted_positive: 64 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.171974786588701 0.173738874441179 0.178919874914129 
##               315                 1                 1                 1 
## 0.196179163902026 0.200093251779467 0.202852777457336 0.210011791381826 
##                 1                 1                 1                 1 
## 0.212974437048088 0.216757518793664 0.220080822230225 0.222802491803171 
##                 1                 1                 1                 1 
## 0.223139990932436 0.223153339029502  0.22342281144256 0.224031896510711 
##                 1                 1                 1                 1 
## 0.224593291408621 0.225327459271123 0.231263133283476 0.232693790570106 
##                 1                 1                 1                 1 
## 0.233230530297115 0.239870080725128 0.239965266489961  0.24016481582862 
##                 1                 1                 1                 1 
## 0.240505104537631 0.242693026068864 0.244919048471134 0.248369771091455 
##                 1                 1                 1                 1 
##   0.2519356440883 0.254280171035428 0.261928931844667 0.265454157532872 
##                 1                 1                 1                 1 
## 0.268616812668484 0.269555193048487 0.271282411941263 0.273649786074466 
##                 1                 1                 1                 1 
## 0.274051701410776 0.276964226221027 0.278035715897949 0.279828285542451 
##                 1                 1                 1                 1 
##  0.28265111043032 0.290884394914459 0.294705364486349 0.298016091682923 
##                 1                 1                 1                 1 
##  0.30344798547644 0.305608509170536  0.30855758972218 0.311533310961849 
##                 1                 1                 1                 1 
## 0.314856703332522 0.320513918748529 0.327073579207257 0.331696198012442 
##                 1                 1                 1                 1 
## 0.331825950371282 0.332169310955188 0.340487914419054 0.343765500005886 
##                 1                 1                 1                 1 
## 0.346296843120759 0.361388157842304  0.38468244237212 0.399717177298722 
##                 1                 1                 1                 1 
## 0.457424954282519   0.4999085540431 0.525316623230191 0.543521048744137 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.250674872011302 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: B.GC.prePB remain_type: B.GC.DZ predicted_positive: 11 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.172529638316355 0.179714547186483 0.181263523571598 
##               369                 1                 1                 1 
## 0.182424928157934 0.183972600002209 0.187610382198632 0.197929921928023 
##                 1                 1                 1                 1 
## 0.223656705938047 0.257156990167656 
##                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.184817849460145 0.241143344115027 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.prePB remain_type: B.GC.LZ predicted_positive: 21 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.153451461580979 0.173875200334031  0.17534617920613 
##               360                 1                 1                 1 
## 0.179301792029225 0.196462171682056  0.19655391994193 0.198878106127906 
##                 1                 1                 1                 1 
## 0.204711490531486 0.208701101190236 0.213940527158371 0.217044072945308 
##                 1                 1                 1                 1 
## 0.225823767963686 0.244330243450618 0.245385824959086 0.248224122468142 
##                 1                 1                 1                 1 
## 0.257447036641531 0.264443254117819 0.322009638036563 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.127699247516295 0.151120191020321 0.276241764931598 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.prePB remain_type: T.CD4..TfH.GC predicted_positive: 4 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.161836382593542 
##               377                 1 
## [1] "non-GC spot"
## 
##                 0 0.158631036956849 0.244082623590589 0.455997827717548 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.GC.prePB remain_type: B.Cycling predicted_positive: 61 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.159453265308961 0.174410579135085 0.176477051676348 
##               318                 1                 1                 1 
## 0.177773979378379 0.182600314269322 0.187533517337576 0.188168395765413 
##                 1                 1                 1                 1 
## 0.189369397938574  0.19517854404028 0.195477070273461 0.202471742007972 
##                 1                 1                 1                 1 
## 0.203720277833395 0.205526814769778 0.208586635289791 0.209622588961652 
##                 1                 1                 1                 1 
## 0.213738632521127 0.215050611837805 0.217284311120625 0.217502536016651 
##                 1                 1                 1                 1 
## 0.220439221907354 0.220701837613259 0.221554608334512 0.226660603538157 
##                 1                 1                 1                 1 
## 0.230523289618702 0.230843836956331 0.231927964990266 0.232928379745833 
##                 1                 1                 1                 1 
## 0.237316639808796 0.237718470203871 0.244672214821238 0.247262820788083 
##                 1                 1                 1                 1 
## 0.256055184107777  0.25715592734591 0.257651722054026 0.259511833675273 
##                 1                 1                 1                 1 
## 0.260681730326315 0.265678028270081 0.266677423414239 0.267281593027704 
##                 1                 1                 1                 1 
## 0.267605031880448 0.273349364385036 0.274334464198432 0.276542572745644 
##                 1                 1                 1                 1 
## 0.282263840974689  0.28377099851252 0.284226011465304 0.293270350138388 
##                 1                 1                 1                 1 
## 0.294085653955301 0.298164506114632 0.315694334967068 0.316551086745508 
##                 1                 1                 1                 1 
##  0.32356479981133 0.323845863733543 0.328519410479785 0.331987716653305 
##                 1                 1                 1                 1 
## 0.345508909548634  0.35178923064267 0.355738946404013 0.375443723456579 
##                 1                 1                 1                 1 
##  0.45578142870729 
##                 1 
## [1] "non-GC spot"
## 
##                 0 0.235741123133085 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: T.CD4..TfH.GC remain_type: B.GC.DZ predicted_positive: 9 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.156685726776613 0.184865488051379   0.1849901244973 
##               371                 1                 1                 1 
## 0.187101330465365 0.225728496847423 0.256053637853188 0.266700561583537 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.273489317696242 0.286330808719648 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: T.CD4..TfH.GC remain_type: B.GC.LZ predicted_positive: 29 actual_positive: 378"
## [1] "GC spot"
## 
##                 0  0.18388204422848 0.187045332177156 0.191228684292098 
##               349                 1                 1                 1 
## 0.195009878950898 0.197118859324834  0.20519272242647 0.212023376550012 
##                 1                 1                 1                 1 
## 0.221108633363886 0.222159960000036 0.225916556461502 0.227210271890659 
##                 1                 1                 1                 1 
## 0.229512798746471 0.242089693286083 0.245725612370704 0.250579679237698 
##                 1                 1                 1                 1 
## 0.257828083489229 0.261946702431095 0.265805903359958 0.273741760411631 
##                 1                 1                 1                 1 
## 0.274706716397681 0.283708460119019 0.284325554231659 0.286311986920238 
##                 1                 1                 1                 1 
## 0.295217417671771 0.302318945344219 0.302847266239228 0.306550931931169 
##                 1                 1                 1                 1 
## 0.375190496215528 0.378085315753311 
##                 1                 1 
## [1] "non-GC spot"
## 
##    0 
## 3657
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: T.CD4..TfH.GC remain_type: B.GC.prePB predicted_positive: 2 actual_positive: 378"
## [1] "GC spot"
## 
##   0 
## 378 
## [1] "non-GC spot"
## 
##                 0  0.12716973215924 0.292022675254996 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: T.CD4..TfH.GC remain_type: B.Cycling predicted_positive: 94 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.190406887899508 0.192772202668369 0.198848253760882 
##               287                 1                 1                 1 
## 0.204299735681061 0.204337876778359 0.204872385759296 0.208924669607153 
##                 1                 1                 1                 1 
## 0.209257053400975 0.211496658257615 0.216031654431688  0.21709762729873 
##                 1                 1                 1                 1 
## 0.220052639105798 0.227497729044572 0.231666812823795 0.231952497429601 
##                 1                 1                 1                 1 
## 0.236248675630032 0.236657507592836 0.237636617419047 0.238234217208579 
##                 1                 1                 1                 1 
## 0.238857994605966 0.238940828282337 0.239897570461943 0.240042364551473 
##                 1                 1                 1                 1 
## 0.240161807870646 0.240762561575044 0.240818117805605 0.243494238239411 
##                 1                 1                 1                 1 
## 0.244079537363697 0.246324557859588 0.249234511797303  0.25084817678529 
##                 1                 1                 1                 1 
## 0.250985955641852 0.252319559526075 0.256696234036554 0.256870385891062 
##                 1                 1                 1                 1 
## 0.257063512137854 0.257687687404347 0.258409767578268  0.25933423170245 
##                 1                 1                 1                 1 
## 0.263932221462069 0.265254662390632 0.266292532950003 0.266449969100418 
##                 1                 1                 1                 1 
## 0.267026153787476 0.268658514303409 0.269388930799184 0.270874524077877 
##                 1                 1                 1                 1 
## 0.272469210417495  0.27250108157968 0.274732115182705 0.276359976979979 
##                 1                 1                 1                 1 
## 0.280739249328736 0.283711217113834 0.287235901614762 0.289727207890419 
##                 1                 1                 1                 1 
## 0.291604640707108 0.294031912031533 0.294595123621729 0.297716249404269 
##                 1                 1                 1                 1 
## 0.298914551493621 0.299119316779642 0.299349426655254 0.299981445520111 
##                 1                 1                 1                 1 
## 0.300615171415523 0.303701153303425 0.308244885587085  0.31524364847699 
##                 1                 1                 1                 1 
##  0.31874801826332 0.319145117466738 0.319880882676789 0.320973481757455 
##                 1                 1                 1                 1 
## 0.327713412839948 0.328084542241143 0.337743251628758  0.33935242043992 
##                 1                 1                 1                 1 
## 0.348072053518412  0.35099580340398 0.352973117563481 0.359001745206167 
##                 1                 1                 1                 1 
## 0.359215798318209 0.365059270562192 0.375311367277843 0.390752217415367 
##                 1                 1                 1                 1 
##  0.39078833994041 0.397818921697722 0.405330162616869 0.418655482248461 
##                 1                 1                 1                 1 
## 0.438241703482384  0.46940673648994 0.481319054635804 0.552136748825949 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.294158560264314 0.344408994447622 0.377983533097753 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: B.Cycling remain_type: B.GC.DZ predicted_positive: 47 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.152466082803455 0.165143110690898 0.167045131589636 
##               335                 1                 1                 1 
## 0.170209340496614 0.170944152037746 0.171715025884391 0.178109258182052 
##                 1                 1                 1                 1 
## 0.181411564876772 0.181947475725192 0.188852123005289 0.191654491509232 
##                 1                 1                 1                 1 
## 0.196221237793657 0.198823108064958 0.199122767659864 0.199741728417562 
##                 1                 1                 1                 1 
## 0.200139851214655 0.202539649884634 0.210268076509364 0.215562516242769 
##                 1                 1                 1                 1 
## 0.220915480674757  0.22154941545643 0.222315628232966 0.226665646736014 
##                 1                 1                 1                 1 
## 0.227766011973137 0.227901385058615 0.229185536067558 0.232429173980613 
##                 1                 1                 1                 1 
## 0.232693634920602 0.235102941370035 0.238209414515897 0.243383955336565 
##                 1                 1                 1                 1 
##  0.24596352628379 0.246088931251677 0.251209744693965 0.257678450452016 
##                 1                 1                 1                 1 
## 0.259923362375512 0.267367731825853 0.280131813825671  0.29235984468347 
##                 1                 1                 1                 1 
## 0.329703128632251 0.332866456522985 0.442186276327341 0.445236386657561 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.165031659936766 0.206099995398573 0.236387838019131 
##              3653                 1                 1                 1 
##  0.29923978433833 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.Cycling remain_type: B.GC.LZ predicted_positive: 195 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.154120191894837 0.154770714269136 0.163246570165205 
##               191                 1                 1                 1 
## 0.171305513245733 0.171727963353476 0.172808029674653 0.174143036655977 
##                 1                 1                 1                 1 
## 0.182259040192259  0.18629501100031 0.187232958093707 0.187617499847599 
##                 1                 1                 1                 1 
## 0.188981162806411 0.190240067257646 0.191733131011865 0.191806435369727 
##                 1                 1                 1                 1 
## 0.192263369711629 0.193859234480602 0.197314502136421 0.198019593208388 
##                 1                 1                 1                 1 
## 0.198062911274047 0.199779103178512 0.203408504874091 0.204353601020329 
##                 1                 1                 1                 1 
## 0.204548757619474 0.205047855260778 0.205640186343609 0.205951827472102 
##                 1                 1                 1                 1 
## 0.206156234494582 0.206431522816466  0.20761618126234 0.208001758968361 
##                 1                 1                 1                 1 
## 0.208121240482366 0.208661446479841 0.208751392108347 0.209581375133838 
##                 1                 1                 1                 1 
## 0.210618709715396 0.210732022642461 0.211061365436784 0.211592813292629 
##                 1                 1                 1                 1 
## 0.211978811572568 0.211986260663022 0.212262233440126 0.212996802760506 
##                 1                 1                 1                 1 
## 0.213064112302004 0.214379364299207 0.214824944974754 0.215121624741083 
##                 1                 1                 1                 1 
## 0.219411156123442 0.220905585314934 0.221028158186776 0.221264180720479 
##                 1                 1                 1                 1 
## 0.222228677804343 0.222638546578683 0.224665973718535  0.22485112434182 
##                 1                 1                 1                 1 
## 0.225154046251023 0.225762792927828 0.226746471153116 0.226970538631381 
##                 1                 1                 1                 1 
## 0.227083829842114 0.227249918785383 0.228194662046463 0.228559355397645 
##                 1                 1                 1                 1 
## 0.229021687451888 0.229386483916639 0.229797744531412 0.230764956565486 
##                 1                 1                 1                 1 
## 0.231167470165428 0.231311099381793 0.231941377645724 0.232023043605396 
##                 1                 1                 1                 1 
## 0.233075001001591 0.233400186175769 0.233813204738548 0.233819976924508 
##                 1                 1                 1                 1 
## 0.234050330234443 0.237321999003186 0.237361443426667  0.23848091308766 
##                 1                 1                 1                 1 
## 0.240169955745592 0.242201903513373 0.242243426974165 0.242833797894394 
##                 1                 1                 1                 1 
## 0.243268959431226 0.245658008281837 0.245787369541636 0.246980750459743 
##                 1                 1                 1                 1 
## 0.247112613662572 0.247860832733844 0.248098671484262 0.248262188965311 
##                 1                 1                 1                 1 
## 0.248911658798774 0.250823293753772 0.251216841159636 0.251628296903398 
##                 1                 1                 1                 1 
## 0.252894853389262 0.253255394879137 0.254273956879187 0.255120597790303 
##                 1                 1                 1                 1 
## 0.256827389178313 0.257299899269001 0.257548277708974 0.258195746490966 
##                 1                 1                 1                 1 
## 0.259472824044361 0.260515095345575 0.261556251670699 0.261719416678369 
##                 1                 1                 1                 1 
## 0.261721090325251 0.264393502937509 0.264709097944007 0.265954510585577 
##                 1                 1                 1                 1 
## 0.266044086530176 0.267461424946679 0.267774325925346 0.268521711565581 
##                 1                 1                 1                 1 
## 0.271112437136655 0.271232413732167 0.271279289295045 0.272300520171792 
##                 1                 1                 1                 1 
## 0.272935391328462 0.273850928552622 0.274284356372332 0.274626645050754 
##                 1                 1                 1                 1 
## 0.276214735372125 0.277159599754589 0.277751954527194 0.278013525383494 
##                 1                 1                 1                 1 
## 0.278243541600636  0.27881828643073 0.281802551864802  0.28183797078674 
##                 1                 1                 1                 1 
## 0.283521572088076 0.284282018235434 0.286945777119376 0.287164271032887 
##                 1                 1                 1                 1 
## 0.287528530936106  0.28773383500615 0.288171399460416 0.288335694692677 
##                 1                 1                 1                 1 
## 0.289257913750241 0.290280308322991 0.291745890395552 0.292312963405264 
##                 1                 1                 1                 1 
##   0.2949544789376 0.296628249364059 0.297401294699599 0.299047569477496 
##                 1                 1                 1                 1 
## 0.299180659648422 0.300910593621139 0.303212694901733 0.304675197158599 
##                 1                 1                 1                 1 
## 0.304831740488945 0.305462190277889 0.306247382376971 0.306538412866291 
##                 1                 1                 1                 1 
## 0.307052403656942 0.309219168318571  0.31099827087515 0.313389976721064 
##                 1                 1                 1                 1 
## 0.318806266465956 0.320952377064593 0.324500956445906 0.326050139589852 
##                 1                 1                 1                 1 
## 0.327663323649924 0.329617580046964 0.329902839145403 0.334639691618395 
##                 1                 1                 1                 1 
## 0.340725365602272 0.341481084531233 0.343054024776343 0.343971296147601 
##                 1                 1                 1                 1 
##  0.34736741902403  0.35917002691358 0.359437236571377 0.360103989669505 
##                 1                 1                 1                 1 
## 0.385028606334825 0.387189127445404  0.39097529454666 0.406313545216238 
##                 1                 1                 1                 1 
## 0.415061084800975 0.415509356253363 0.425853972526709 0.428798561803184 
##                 1                 1                 1                 1 
##  0.44983613272734 0.472793159944533 0.502969131714749 0.539567895990498 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.168633420257567 0.205037843635469 0.211258957584836 
##              3649                 1                 1                 1 
## 0.267556087291319 0.286603116851398  0.30896419130652 0.309340285482271 
##                 1                 1                 1                 1 
##  0.33264034783165 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.Cycling remain_type: B.GC.prePB predicted_positive: 3 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.147581690179008 
##               377                 1 
## [1] "non-GC spot"
## 
##                 0 0.106830135266293 0.184298933242049 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: B.Cycling remain_type: T.CD4..TfH.GC predicted_positive: 5 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.186262796047352 0.191359640733918 0.291131121068923 
##               375                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0  0.23995847417867 0.428162821436838 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: FDC remain_type: B.GC.DZ predicted_positive: 22 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.188615670287114 0.203098968232173 0.216539601186953 
##               360                 1                 1                 1 
## 0.242634703610394 0.247565073396384 0.252217544762809 0.262654442620813 
##                 1                 1                 1                 1 
##  0.27289182043013 0.277201472814569 0.286637957525055 0.288150016664362 
##                 1                 1                 1                 1 
## 0.290997111918552 0.292044190136111 0.293516764210672 0.297573861144895 
##                 1                 1                 1                 1 
## 0.322003425092026 0.325137474677073 0.337337826401395 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.230111715150295 0.263081334808753 0.282892976856037 
##              3653                 1                 1                 1 
## 0.305683013880545 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: FDC remain_type: B.GC.LZ predicted_positive: 85 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.172872651025256 0.223934206111431 0.235709601302832 
##               296                 1                 1                 1 
## 0.243934788607791 0.245552772311699 0.255013799400214 0.260139988475658 
##                 1                 1                 1                 1 
## 0.261186444098397 0.262014120344372 0.267485098722252 0.273842727713542 
##                 1                 1                 1                 1 
## 0.280975256154964 0.284602258588947 0.293736468434734 0.295412100528589 
##                 1                 1                 1                 1 
## 0.295681495690994 0.296499618400245 0.298645160824064 0.300007084528795 
##                 1                 1                 1                 1 
## 0.301702925757404 0.307780375829258 0.312078765158544 0.315371518498901 
##                 1                 1                 1                 1 
## 0.315587465676779 0.319102444152837 0.321078482479879 0.322961172691596 
##                 1                 1                 1                 1 
## 0.323953205207702 0.324243789303799 0.326033697150121 0.327841054724984 
##                 1                 1                 1                 1 
## 0.329133601337792 0.330934749155413  0.33191359776576 0.341161146481493 
##                 1                 1                 1                 1 
## 0.344284885683902 0.350896447145397 0.351758488272907 0.352054588314258 
##                 1                 1                 1                 1 
## 0.354248194524486 0.355392095729372 0.357279117431815 0.357860675867576 
##                 1                 1                 1                 1 
## 0.358481186734607 0.361370777763656 0.361542480568916 0.362664976316538 
##                 1                 1                 1                 1 
## 0.373025112560261 0.374726159510252 0.378411313551809  0.38188567737563 
##                 1                 1                 1                 1 
## 0.384071433681411 0.387424651600797  0.38918186871385 0.393291558862523 
##                 1                 1                 1                 1 
## 0.393424856997987 0.394266622978195 0.394912637551516 0.395826926243019 
##                 1                 1                 1                 1 
## 0.401460261882615 0.402073940323462 0.403053143306733 0.405233958861828 
##                 1                 1                 1                 1 
## 0.405381871646693  0.40879892760054 0.429086309540976 0.434639150594452 
##                 1                 1                 1                 1 
## 0.439161662333576 0.448418130577419 0.450571734593622 0.451118918185302 
##                 1                 1                 1                 1 
## 0.456653837433308 0.457329293068952 0.461735789904595 0.462471574501915 
##                 1                 1                 1                 1 
## 0.466272254533456 0.475064750156287 0.478450667706394 0.484287199652783 
##                 1                 1                 1                 1 
## 0.515749747610999 0.578981823959023 0.591059663357023 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.152109358043293 0.179282346273112 0.309410732493002 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: FDC remain_type: B.GC.prePB predicted_positive: 3 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.241951112844146 
##               377                 1 
## [1] "non-GC spot"
## 
##                 0 0.201394450052036                 1 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: FDC remain_type: T.CD4..TfH.GC predicted_positive: 8 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.279056658645284 0.284562719673132 0.370296501073489 
##               374                 1                 1                 1 
## 0.460621970456405 
##                 1 
## [1] "non-GC spot"
## 
##                 0 0.191716582374861 0.296758486110963 0.316801299962947 
##              3653                 1                 1                 1 
## 0.460516222910319 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: FDC remain_type: B.Cycling predicted_positive: 106 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.228199613972335 0.229771775835529 0.239080000317213 
##               276                 1                 1                 1 
## 0.239518128206818  0.26066785673831 0.268648354976798 0.268895777199184 
##                 1                 1                 1                 1 
## 0.270378522596882 0.270496530595934 0.279497700146883 0.281364764885497 
##                 1                 1                 1                 1 
## 0.282104051308028 0.286750786802994 0.287516447470597 0.289517067848898 
##                 1                 1                 1                 1 
##  0.29173968642236 0.292242849341343 0.293258847879852 0.294653951213074 
##                 1                 1                 1                 1 
## 0.294760061887752 0.296732704873549 0.298861474616011 0.304388520204193 
##                 1                 1                 1                 1 
## 0.306924214242904 0.307640159270768 0.308289212739275 0.311446848984689 
##                 1                 1                 1                 1 
## 0.312826894900991 0.312907430840096 0.314959618583764 0.323462724597811 
##                 1                 1                 1                 1 
## 0.330867147759036 0.331698512638679 0.332867526699793 0.333194095977528 
##                 1                 1                 1                 1 
## 0.338266254836451 0.342034246916293 0.343381747589153 0.349001143175949 
##                 1                 1                 1                 1 
## 0.349079185799141 0.359753525925345 0.360897289102133 0.366472925004649 
##                 1                 1                 1                 1 
## 0.366882745498839 0.367661255324435 0.367663350289474 0.368808081489969 
##                 1                 1                 1                 1 
## 0.369127720264939 0.370328260981377 0.372981729650895 0.381229961165465 
##                 1                 1                 1                 1 
## 0.382570943697421 0.383107615052256 0.383546251847073 0.385354698128649 
##                 1                 1                 1                 1 
## 0.388970826332177 0.391706198504696 0.397995563102324 0.399724366892045 
##                 1                 1                 1                 1 
##  0.40315385605693  0.40325682649599 0.403589108515524 0.404351315864629 
##                 1                 1                 1                 1 
##  0.40754461896285 0.413075106013046 0.413285466257068 0.414348213255251 
##                 1                 1                 1                 1 
## 0.417769090515353 0.419411693766385 0.422988715774875 0.424527201973641 
##                 1                 1                 1                 1 
## 0.425983748229367 0.430413574510555 0.432256685101539  0.43631497031131 
##                 1                 1                 1                 1 
## 0.436861465896173 0.438365074593513 0.448129710395983 0.448235025639166 
##                 1                 1                 1                 1 
## 0.449808621415865  0.45103978382777 0.451384910206497 0.452497844123753 
##                 1                 1                 1                 1 
## 0.454456692440559 0.459187012074392 0.461979569152812 0.464335606126144 
##                 1                 1                 1                 1 
## 0.464745294829317 0.469258459035845 0.471744145133214 0.475077879486644 
##                 1                 1                 1                 1 
## 0.477183669989285 0.485625933435227 0.505917384618904 0.506503363020975 
##                 1                 1                 1                 1 
## 0.509952468293514 0.511719538229478 0.512747214733329 0.514214186917832 
##                 1                 1                 1                 1 
## 0.524874700812831 0.544959247947232 0.555529650393032 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.387402972167572 0.399452366617143  0.41533873383275 
##              3653                 1                 1                 1 
## 0.519735652804261 
##                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 5 rows containing missing values (geom_point).

## [1] "deleted_type: complete remain_type: B.GC.DZ predicted_positive: 8 actual_positive: 378"
## [1] "GC spot"
## 
##                 0  0.18765638605767 0.196891700357968 0.211171845142658 
##               372                 1                 1                 1 
## 0.220450551695887 0.222234071509082 0.319294868592376 
##                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0  0.19220386159362 0.223536133718597 
##              3655                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: complete remain_type: B.GC.LZ predicted_positive: 29 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.180354269397146 0.180676445245764 0.184065343421816 
##               350                 1                 1                 1 
## 0.190665917080471 0.196155297183303 0.203944777009589 0.205800671046546 
##                 1                 1                 1                 1 
## 0.206700704575008 0.211867200167801 0.218754855851122 0.230263338649266 
##                 1                 1                 1                 1 
## 0.232879767232675 0.240188727509414 0.240514978420727  0.24771207426435 
##                 1                 1                 1                 1 
## 0.249534796481811 0.253056651206707 0.262876704329313 0.266087189542783 
##                 1                 1                 1                 1 
## 0.275200177163388 0.281751395143975 0.284393066079981 0.289396923431834 
##                 1                 1                 1                 1 
## 0.302557720637076 0.309370119052651 0.395066489661783 0.480422463222624 
##                 1                 1                 1                 1 
## 0.529062286012755 
##                 1 
## [1] "non-GC spot"
## 
##                 0 0.235088503200021 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: complete remain_type: B.GC.prePB predicted_positive: 1 actual_positive: 378"
## [1] "GC spot"
## 
##   0 
## 378 
## [1] "non-GC spot"
## 
##                 0 0.151514267571877 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: complete remain_type: T.CD4..TfH.GC predicted_positive: 5 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.186292663238993 0.244343731180776 
##               376                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.183334880259061 0.196820352480787 0.385968989303328 
##              3654                 1                 1                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## [1] "deleted_type: complete remain_type: B.Cycling predicted_positive: 52 actual_positive: 378"
## [1] "GC spot"
## 
##                 0 0.195364023321049  0.20438517229301 0.205842064320721 
##               327                 1                 1                 1 
## 0.205923848787063  0.22030637703886 0.221994935535965  0.22377293384192 
##                 1                 1                 1                 1 
## 0.224798168916443 0.227161092897616 0.228869129443005 0.229585396174736 
##                 1                 1                 1                 1 
## 0.235236449559865  0.23765131206781 0.239251082014736 0.242172837602361 
##                 1                 1                 1                 1 
##  0.24302033034494 0.249351822688749 0.259763402769832 0.262571014587541 
##                 1                 1                 1                 1 
##  0.27230974034403 0.275940808854295 0.277031453220096 0.284032545468245 
##                 1                 1                 1                 1 
## 0.285297553026466 0.285805463948311 0.287493365172948 0.288780950453101 
##                 1                 1                 1                 1 
## 0.295710009482336 0.298709749272331 0.300289975511387 0.301411641141957 
##                 1                 1                 1                 1 
## 0.302951601437253  0.30437341376777  0.30808210776181 0.312126811072563 
##                 1                 1                 1                 1 
## 0.313977541093485 0.331872481710262 0.334216258312612 0.335308244367928 
##                 1                 1                 1                 1 
## 0.335846462380555 0.338278439569191 0.346344857387304 0.348424186414556 
##                 1                 1                 1                 1 
## 0.354544466565788 0.384162275911473 0.389336206246645 0.465803439863011 
##                 1                 1                 1                 1 
## 0.498375730330153 0.500561356503688 0.519577536777376                 1 
##                 1                 1                 1                 1 
## [1] "non-GC spot"
## 
##                 0 0.443140605304528 
##              3656                 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 6 rows containing missing values (geom_point).

df <- data.frame(deleted_cell_type=names(average_PR_curve.area),
                average_area_PR_curve=average_PR_curve.area, similarity = "similar"  )
df[df$deleted_cell_type %in% c("B.Cycling" , "FDC"),"similarity"] = "distinct"

ref_line = df["complete",2]
df = df[df$deleted_cell_type!="complete",]

p <-ggplot(data=df, aes(x=deleted_cell_type,y=average_area_PR_curve,fill=similarity)) +  geom_bar(stat="identity") + geom_hline(yintercept=ref_line, linetype="dashed", color = "red")
p

Session information

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ROCR_1.0-11           SPOTlight_0.99.4      RColorBrewer_1.1-2   
##  [4] tidyr_1.2.0           tibble_3.1.6          stringr_1.4.0        
##  [7] bigmemory_4.5.36      Biobase_2.52.0        BiocGenerics_0.38.0  
## [10] ggplot2_3.3.5         dplyr_1.0.8           SeuratDisk_0.0.0.9019
## [13] SeuratObject_4.0.4    Seurat_4.1.0          data.table_1.14.2    
## 
## loaded via a namespace (and not attached):
##   [1] NMF_0.23.0                  plyr_1.8.6                 
##   [3] igraph_1.2.11               lazyeval_0.2.2             
##   [5] splines_4.1.0               listenv_0.8.0              
##   [7] scattermore_0.7             GenomeInfoDb_1.28.4        
##   [9] gridBase_0.4-7              digest_0.6.29              
##  [11] foreach_1.5.2               htmltools_0.5.2            
##  [13] fansi_1.0.2                 magrittr_2.0.2             
##  [15] tensor_1.5                  cluster_2.1.2              
##  [17] doParallel_1.0.17           globals_0.14.0             
##  [19] matrixStats_0.61.0          spatstat.sparse_2.1-0      
##  [21] colorspace_2.0-2            ggrepel_0.9.1              
##  [23] xfun_0.29                   RCurl_1.98-1.6             
##  [25] crayon_1.4.2                jsonlite_1.7.3             
##  [27] bigmemory.sri_0.1.3         scatterpie_0.1.7           
##  [29] spatstat.data_2.1-2         survival_3.2-13            
##  [31] zoo_1.8-9                   iterators_1.0.14           
##  [33] glue_1.6.1                  polyclip_1.10-0            
##  [35] registry_0.5-1              gtable_0.3.0               
##  [37] nnls_1.4                    zlibbioc_1.38.0            
##  [39] XVector_0.32.0              leiden_0.3.9               
##  [41] DelayedArray_0.18.0         SingleCellExperiment_1.14.1
##  [43] future.apply_1.8.1          abind_1.4-5                
##  [45] scales_1.1.1                DBI_1.1.2                  
##  [47] rngtools_1.5.2              miniUI_0.1.1.1             
##  [49] Rcpp_1.0.8                  viridisLite_0.4.0          
##  [51] xtable_1.8-4                reticulate_1.24            
##  [53] spatstat.core_2.3-2         bit_4.0.4                  
##  [55] stats4_4.1.0                htmlwidgets_1.5.4          
##  [57] httr_1.4.2                  ellipsis_0.3.2             
##  [59] ica_1.0-2                   pkgconfig_2.0.3            
##  [61] farver_2.1.0                sass_0.4.0                 
##  [63] uwot_0.1.11                 deldir_1.0-6               
##  [65] utf8_1.2.2                  tidyselect_1.1.1           
##  [67] labeling_0.4.2              rlang_1.0.1                
##  [69] reshape2_1.4.4              later_1.3.0                
##  [71] munsell_0.5.0               tools_4.1.0                
##  [73] cli_3.1.1                   generics_0.1.2             
##  [75] ggridges_0.5.3              evaluate_0.14              
##  [77] fastmap_1.1.0               yaml_2.2.2                 
##  [79] goftest_1.2-3               knitr_1.37                 
##  [81] bit64_4.0.5                 fitdistrplus_1.1-6         
##  [83] purrr_0.3.4                 RANN_2.6.1                 
##  [85] pbapply_1.5-0               future_1.23.0              
##  [87] nlme_3.1-152                mime_0.12                  
##  [89] hdf5r_1.3.5                 compiler_4.1.0             
##  [91] rstudioapi_0.13             plotly_4.10.0              
##  [93] png_0.1-7                   spatstat.utils_2.3-0       
##  [95] tweenr_1.0.2                bslib_0.3.1                
##  [97] stringi_1.7.6               highr_0.9                  
##  [99] lattice_0.20-44             Matrix_1.4-0               
## [101] vctrs_0.3.8                 pillar_1.7.0               
## [103] lifecycle_1.0.1             spatstat.geom_2.3-1        
## [105] lmtest_0.9-39               jquerylib_0.1.4            
## [107] RcppAnnoy_0.0.19            bitops_1.0-7               
## [109] cowplot_1.1.1               irlba_2.3.5                
## [111] GenomicRanges_1.44.0        httpuv_1.6.5               
## [113] patchwork_1.1.1             R6_2.5.1                   
## [115] promises_1.2.0.1            KernSmooth_2.23-20         
## [117] gridExtra_2.3               IRanges_2.26.0             
## [119] parallelly_1.30.0           codetools_0.2-18           
## [121] MASS_7.3-54                 assertthat_0.2.1           
## [123] SummarizedExperiment_1.22.0 pkgmaker_0.32.2            
## [125] withr_2.4.3                 sctransform_0.3.3          
## [127] GenomeInfoDbData_1.2.6      S4Vectors_0.30.2           
## [129] mgcv_1.8-38                 ggfun_0.0.5                
## [131] grid_4.1.0                  rpart_4.1-15               
## [133] rmarkdown_2.11              MatrixGenerics_1.4.3       
## [135] Rtsne_0.15                  ggforce_0.3.3              
## [137] shiny_1.7.1